import numpy as np
import pandas as pd
from scipy import stats
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt
import seaborn as sns

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 读取Excel数据
df = pd.read_excel('含达标时间的男胎数据.xlsx')

# 提取需要的列并重命名
df = df[['孕妇BMI（处理后）', '达标时间']]
df.columns = ['BMI', '达标时间']

# 定义BMI分组
bmi_groups = [
    {'min': 26.619, 'max': 29.891, 'name': '组1', 'mean_bmi': 28.926, 'n': 59},
    {'min': 29.969, 'max': 32.136, 'name': '组2', 'mean_bmi': 30.933, 'n': 81},
    {'min': 32.311, 'max': 35.059, 'name': '组3', 'mean_bmi': 33.404, 'n': 66},
    {'min': 35.380, 'max': 39.302, 'name': '组4', 'mean_bmi': 36.811, 'n': 28}
]

# 医疗风险函数
def risk_function(t):
    if t <= 12:
        return 0.1
    elif t <= 27:
        return 0.1 + 0.03 * (t - 12)
    else:
        return 1.0

# 使用核密度估计拟合达标时间分布
def fit_kde(data, bw_method='silverman'):
    kde = stats.gaussian_kde(data, bw_method=bw_method)
    return kde

# 计算达标比例函数
def cumulative_proportion(kde, t):
    # 计算从10周到t周的累积概率
    x = np.linspace(10, t, 1000)
    y = kde(x)
    return np.trapz(y, x)

# 目标函数：最大化 P_d(t) * (1 - R(t))
def objective_function(t, kde):
    P_d = cumulative_proportion(kde, t)
    R_t = risk_function(t)
    return -P_d * (1 - R_t)  # 负号是因为我们要最大化

# 为每个分组计算最佳检测时点
results = []
group_stats = []  # 存储每组的统计信息

for group in bmi_groups:
    # 筛选该组的数据
    group_data = df[(df['BMI'] >= group['min']) & (df['BMI'] <= group['max'])]
    
    # 检查是否有足够的数据
    if len(group_data) < 5:
        print(f"警告: {group['name']} 的数据量不足 ({len(group_data)}个样本)，跳过计算")
        continue
        
    达标时间数据 = group_data['达标时间'].values
    
    # 计算描述性统计
    stats_dict = {
        '分组': group['name'],
        'BMI范围': f"{group['min']:.1f}-{group['max']:.1f}",
        '平均BMI': group['mean_bmi'],
        '样本数': len(group_data),
        '达标时间均值': np.mean(达标时间数据),
        '达标时间中位数': np.median(达标时间数据),
        '达标时间标准差': np.std(达标时间数据),
        '达标时间最小值': np.min(达标时间数据),
        '达标时间最大值': np.max(达标时间数据),
        '25%分位数': np.percentile(达标时间数据, 25),
        '75%分位数': np.percentile(达标时间数据, 75)
    }
    group_stats.append(stats_dict)
    
    # 拟合核密度估计
    try:
        kde = fit_kde(达标时间数据)
    except Exception as e:
        print(f"错误: 无法为 {group['name']} 拟合核密度估计: {e}")
        continue
    
    # 优化目标函数
    try:
        res = minimize_scalar(
            objective_function, 
            args=(kde,),
            bounds=(10, 25),
            method='bounded'
        )
        
        # 提取结果
        optimal_t = res.x
        P_d = cumulative_proportion(kde, optimal_t)
        R_t = risk_function(optimal_t)
        objective_value = -res.fun  # 恢复正值
        
        # 计算几个关键时间点的值
        t_16 = 16
        P_d_16 = cumulative_proportion(kde, t_16)
        R_t_16 = risk_function(t_16)
        objective_16 = P_d_16 * (1 - R_t_16)
        
        t_17 = 17
        P_d_17 = cumulative_proportion(kde, t_17)
        R_t_17 = risk_function(t_17)
        objective_17 = P_d_17 * (1 - R_t_17)
        
        t_18 = 18
        P_d_18 = cumulative_proportion(kde, t_18)
        R_t_18 = risk_function(t_18)
        objective_18 = P_d_18 * (1 - R_t_18)
        
        results.append({
            '分组': group['name'],
            'BMI范围': f"[{group['min']:.3f}, {group['max']:.3f}]",
            '样本数': len(group_data),
            '最佳检测时点(周)': f"{optimal_t:.1f}",
            '预期达标比例': f"{P_d*100:.1f}%",
            '风险系数': f"{R_t:.2f}",
            '目标函数值': f"{objective_value:.3f}",
            '16周目标值': f"{objective_16:.3f}",
            '17周目标值': f"{objective_17:.3f}",
            '18周目标值': f"{objective_18:.3f}"
        })
    except Exception as e:
        print(f"错误: 无法为 {group['name']} 优化目标函数: {e}")

# 显示结果
if results:
    results_df = pd.DataFrame(results)
    print("最佳NIPT检测时点优化结果:")
    print(results_df.to_string(index=False))
    
    # 显示描述性统计
    stats_df = pd.DataFrame(group_stats)
    print("\n各分组达标时间的描述性统计:")
    print(stats_df.to_string(index=False))
    
    # 可视化结果
    fig, axes = plt.subplots(2, 2, figsize=(15, 12))
    
    # 1. 目标函数曲线
    for i, group in enumerate(bmi_groups):
        group_data = df[(df['BMI'] >= group['min']) & (df['BMI'] <= group['max'])]
        
        if len(group_data) < 5:
            continue
            
        达标时间数据 = group_data['达标时间'].values
        kde = fit_kde(达标时间数据)
        
        t_values = np.linspace(10, 25, 100)
        objective_values = [-objective_function(t, kde) for t in t_values]
        
        axes[0, 0].plot(t_values, objective_values, label=f"{group['name']} (BMI: {group['mean_bmi']:.1f})")
        
        optimal_t = float([r for r in results if r['分组'] == group['name']][0]['最佳检测时点(周)'])
        optimal_value = float([r for r in results if r['分组'] == group['name']][0]['目标函数值'])
        axes[0, 0].plot(optimal_t, optimal_value, 'ro')
        axes[0, 0].text(optimal_t, optimal_value, f"{optimal_t:.1f}周", ha='left', va='bottom')

    axes[0, 0].set_xlabel('检测孕周')
    axes[0, 0].set_ylabel('目标函数值 (P_d(t) * (1 - R(t)))')
    axes[0, 0].set_title('各BMI分组的最佳NIPT检测时点')
    axes[0, 0].legend()
    axes[0, 0].grid(True)
    
    # 2. 达标时间分布直方图和核密度估计
    colors = ['blue', 'green', 'red', 'orange']
    
    for i, group in enumerate(bmi_groups):
        group_data = df[(df['BMI'] >= group['min']) & (df['BMI'] <= group['max'])]
        
        if len(group_data) < 5:
            continue
            
        # 绘制直方图
        axes[0, 1].hist(group_data['达标时间'], bins=20, alpha=0.3, 
                label=f"{group['name']} (n={len(group_data)})", color=colors[i])
        
        # 绘制核密度估计曲线
        达标时间数据 = group_data['达标时间'].values
        kde = fit_kde(达标时间数据)
        x = np.linspace(达标时间数据.min(), 达标时间数据.max(), 1000)
        y = kde(x)
        axes[0, 1].plot(x, y * len(group_data) * 1.5, color=colors[i], linewidth=2)
    
    axes[0, 1].set_xlabel('达标时间(周)')
    axes[0, 1].set_ylabel('频数')
    axes[0, 1].set_title('各BMI分组的达标时间分布及核密度估计')
    axes[0, 1].legend()
    axes[0, 1].grid(True)
    
    # 3. 风险函数和达标比例曲线
    t_values = np.linspace(10, 25, 100)
    risk_values = [risk_function(t) for t in t_values]
    axes[1, 0].plot(t_values, risk_values, 'r-', label='风险函数 R(t)', linewidth=2)
    
    for i, group in enumerate(bmi_groups):
        group_data = df[(df['BMI'] >= group['min']) & (df['BMI'] <= group['max'])]
        
        if len(group_data) < 5:
            continue
            
        达标时间数据 = group_data['达标时间'].values
        kde = fit_kde(达标时间数据)
        
        cumulative_values = [cumulative_proportion(kde, t) for t in t_values]
        axes[1, 0].plot(t_values, cumulative_values, label=f"{group['name']} 达标比例", linewidth=2)
    
    axes[1, 0].set_xlabel('孕周')
    axes[1, 0].set_ylabel('值')
    axes[1, 0].set_title('风险函数和各组达标比例随孕周变化')
    axes[1, 0].legend()
    axes[1, 0].grid(True)
    
    # 4. BMI与达标时间的散点图
    axes[1, 1].scatter(df['BMI'], df['达标时间'], alpha=0.5)
    
    for group in bmi_groups:
        axes[1, 1].axvline(x=group['min'], color='r', linestyle='--', alpha=0.5)
        axes[1, 1].axvline(x=group['max'], color='r', linestyle='--', alpha=0.5)
        axes[1, 1].text((group['min']+group['max'])/2, 25, group['name'], 
                horizontalalignment='center', fontsize=10)
    
    axes[1, 1].set_xlabel('BMI')
    axes[1, 1].set_ylabel('达标时间(周)')
    axes[1, 1].set_title('BMI与达标时间的关系')
    axes[1, 1].grid(True)
    
    plt.tight_layout()
    plt.savefig('综合分析.png', dpi=300)
    plt.show()
    
    # 特别关注组1和组2的比较
    print("\n组1和组2的详细比较:")
    group1_data = df[(df['BMI'] >= bmi_groups[0]['min']) & (df['BMI'] <= bmi_groups[0]['max'])]
    group2_data = df[(df['BMI'] >= bmi_groups[1]['min']) & (df['BMI'] <= bmi_groups[1]['max'])]
    
    print(f"组1达标时间均值: {np.mean(group1_data['达标时间']):.2f}周")
    print(f"组2达标时间均值: {np.mean(group2_data['达标时间']):.2f}周")
    
    # 进行t检验
    t_stat, p_value = stats.ttest_ind(group1_data['达标时间'], group2_data['达标时间'])
    print(f"t检验结果: t={t_stat:.3f}, p={p_value:.3f}")
    
    if p_value < 0.05:
        print("两组达标时间有显著差异")
    else:
        print("两组达标时间无显著差异")
    
    # 检查数据分布形状
    print(f"\n组1达标时间偏度: {stats.skew(group1_data['达标时间']):.3f}")
    print(f"组2达标时间偏度: {stats.skew(group2_data['达标时间']):.3f}")
    
    print(f"组1达标时间峰度: {stats.kurtosis(group1_data['达标时间']):.3f}")
    print(f"组2达标时间峰度: {stats.kurtosis(group2_data['达标时间']):.3f}")
    
else:
    print("没有成功计算任何分组的结果")
